clear; close all; clc;

%% load results

cdir = pwd;

idcase = 2; % idcase = 1 for Argentina, idcase = 2 for Spain

savef = 1;               % savef = 1 to save plots
fdir  = 'Figures/IRF';  % location to save plots


if idcase == 1
    cname = 'Argentina';
elseif idcase == 2
    cname = 'Spain';
else
    disp('set idcase=1 for Argenta or idcase=2 for Spain')
    messi
end
    

load([cname,'/OUTPUT_1/prms.txt'])

%--grids
load([cname,'/OUTPUT_1/bvec.txt'])
load([cname,'/OUTPUT_1/zvec.txt'])
load([cname,'/OUTPUT_1/Pzvec.txt'])

%--policies
load([cname,'/OUTPUT_1/bpol.txt'])
load([cname,'/OUTPUT_1/dpol.txt'])
load([cname,'/OUTPUT_1/bopt_loc.txt'])
load([cname,'/OUTPUT_1/Rsched.txt'])

%prms = [Rs,dlta,ybarL,ybarH,kppa,pgL,pgH,bta,gL,gH,sdz,psB]
Rs  = prms(1) ; dlta = prms(2) ; ybarL = prms(3) ; ybarH = prms(4) ; kppa  = prms(5);
pgL = prms(6) ; pgH  = prms(7) ; bta   = prms(8) ; gL    = prms(9) ; gH    = prms(10); 
sdz = prms(11); psB   = prms(12);

psG   = 1.0-psB;

gvec = [gL, gH];
Pg   = [pgL  , 1-pgL; ...
        1-pgH, pgH];
Ps   = [psB  , 1-psB; ...
        1-psG, psG];

Nb = numel(bvec); Nz = numel(zvec); Ng = 2; Ns = 2;

%% computation and reshape

bpolx = bpol; dpolx = dpol; bopt_locx = bopt_loc; Rschedx = Rsched;

bpol     = 189*ones(Nb,Ng,Ns,Nz);
dpol     = 189*ones(Nb,Ng,Ns,Nz);
bopt_loc = 189*ones(Nb,Ng,Ns,Nz);

Rsched = 189*ones(Nb,Ng,Ns);
nsched = 189*ones(Nb,Nb,Ng,Ns);

%--reshape
inn = 0;
for ig = 1:Ng
for is = 1:Ns
for iz = 1:Nz
    for ib = 1:Nb
        inn = inn+1;
        bpol(ib,ig,is,iz)     = bpolx(inn);
        dpol(ib,ig,is,iz)     = dpolx(inn);
        bopt_loc(ib,ig,is,iz) = bopt_locx(inn);
    end
end
end
end

inn = 0;
for ig = 1:Ng
for is = 1:Ns
    for ibnp = 1:Nb
        inn = inn+1;
        Rsched(ibnp,ig,is) = Rschedx(inn);
    end
end
end
for ig = 1:Ng
for is = 1:Ns
    for ib = 1:Nb
        gg = gvec(ig);
        nsched(:,ib,ig,is) = (gg*bvec - (1-dlta)*bvec(ib))./(dlta*Rsched(:,ig,is));
    end
end
end
rhosched = dlta*Rsched-dlta;
rhosched = (rhosched-(Rs-1));

nopt = 189*ones(Nb,Ng,Ns,Nz); nplot = 189*ones(Nb,Ng,Ns,Nz);
ropt = 189*ones(Nb,Ng,Ns,Nz); rplot = 189*ones(Nb,Ng,Ns,Nz);
for ib = 1:Nb
for ig = 1:Ng
for is = 1:Ns
for iz = 1:Nz    
    ibnp = bopt_loc(ib,ig,is,iz); dpx = dpol(ib,ig,is,iz);
    nopt(ib,ig,is,iz) = nsched(ibnp,ib,ig,is);  
    ropt(ib,ig,is,iz) = rhosched(ibnp,ig,is);   
    if dpx == 0
        nplot(ib,ig,is,iz) = nopt(ib,ig,is,iz);
        rplot(ib,ig,is,iz) = ropt(ib,ig,is,iz);
    else
        nplot(ib,ig,is,iz) = NaN;
        rplot(ib,ig,is,iz) = NaN;
    end
end
end
end
end



%% IRF simulation

Tsim = 4; Nsim = 50000;  cix = 1.00;

%--------------------------------------------------------
%---CDF for draws
%-z dist
CDFz = 189*ones(Nz,1);
CDFz(1) = Pzvec(1);
for iz = 2:Nz
    CDFz(iz) = CDFz(iz-1) + Pzvec(iz);
end

%-g dist
CDFg = 189*ones(Ng,Ng);
for iglp = 1:Ng
    CDFg(iglp,1) = Pg(iglp,1);
    for igtp = 2:Ng
        CDFg(iglp,igtp) = CDFg(iglp,igtp-1) + Pg(iglp,igtp);
    end
end

%-s dist
CDFs = 189*ones(Ns,Ns);
for islp = 1:Ns
    CDFs(islp,1) = Ps(islp,1);
    for istp = 2:Ns
        CDFs(islp,istp) = CDFs(islp,istp-1) + Ps(islp,istp);
    end
end
%--------------------------------------------------------


%--------------------------------------------------------
%---initial point
if idcase == 1  % ARG case
    ib0 = find(bvec<0.20,1,'last'); bval  = '20';
elseif  idcase == 2  % SPA case
    ib0 = find(bvec<0.15,1,'last'); bval  = '15';
end
ig0 = 1; iz0 = 9; %floor(Nz/2); is0 = 1;

fsim0 = 0.50;

nn = squeeze(nopt(ib0,ig0,:,iz0))';
rr = squeeze(ropt(ib0,ig0,:,iz0))';
dd = squeeze(dpol(ib0,ig0,:,iz0))';
%--------------------------------------------------------

%--------------------------------------------------------
%---simulaiton: case s_1 = s_B
is0 = 1;
run run_sim.m

IRF_rr_sB = IRF_rr; IRF_rr_sd_sB = IRF_rr_sd; IRF_rr_lb_sB = IRF_rr_sB-cix*IRF_rr_sd_sB; IRF_rr_ub_sB = IRF_rr_sB+cix*IRF_rr_sd_sB;
IRF_nn_sB = IRF_nn; IRF_nn_sd_sB = IRF_nn_sd; IRF_nn_lb_sB = IRF_nn_sB-cix*IRF_nn_sd_sB; IRF_nn_ub_sB = IRF_nn_sB+cix*IRF_nn_sd_sB;
IRF_bb_sB = IRF_bb; IRF_bb_sd_sB = IRF_bb_sd; IRF_bb_lb_sB = IRF_bb_sB-cix*IRF_bb_sd_sB; IRF_bb_ub_sB = IRF_bb_sB+cix*IRF_bb_sd_sB;
IRF_ff_sB = IRF_ff; IRF_ff_sd_sB = IRF_ff_sd; IRF_ff_lb_sB = IRF_ff_sB-cix*IRF_ff_sd_sB; IRF_ff_ub_sB = IRF_ff_sB+cix*IRF_ff_sd_sB;
nuseB = sum(dsim(:,Tsim)==0);
disp(['sims sB with no default = ',num2str(nuseB)])
%--------------------------------------------------------

%--------------------------------------------------------
%---simulaiton: case s_1 = s_G
is0 = 2;

run run_sim.m

IRF_rr_sG = IRF_rr; IRF_rr_sd_sG = IRF_rr_sd; IRF_rr_lb_sG = IRF_rr_sG-cix*IRF_rr_sd_sG; IRF_rr_ub_sG = IRF_rr_sG+cix*IRF_rr_sd_sG;
IRF_nn_sG = IRF_nn; IRF_nn_sd_sG = IRF_nn_sd; IRF_nn_lb_sG = IRF_nn_sG-cix*IRF_nn_sd_sG; IRF_nn_ub_sG = IRF_nn_sG+cix*IRF_nn_sd_sG;
IRF_bb_sG = IRF_bb; IRF_bb_sd_sG = IRF_bb_sd; IRF_bb_lb_sG = IRF_bb_sG-cix*IRF_bb_sd_sG; IRF_bb_ub_sG = IRF_bb_sG+cix*IRF_bb_sd_sG;
IRF_ff_sG = IRF_ff; IRF_ff_sd_sG = IRF_ff_sd; IRF_ff_lb_sG = IRF_ff_sG-cix*IRF_ff_sd_sG; IRF_ff_ub_sG = IRF_ff_sG+cix*IRF_ff_sd_sG;
nuseG = sum(dsim(:,Tsim)==0);
disp(['sims sG with no default = ',num2str(nuseG)])
%--------------------------------------------------------

%% simulation plot

linesz   = 0.5;
labelsz  = 9;
legendsz = 9;
numsz    = 9;
textsz   = 9;

tvec = (1:Tsim)'-1;


fig = figure(202201); clf
plot(tvec,100*IRF_rr_sB,'color',[0.8 0.4 0.4],'LineWidth',linesz,'LineStyle','--') ; hold on
plot(tvec,100*IRF_rr_sG,'color',[0.8 0.4 0.4],'LineWidth',linesz,'LineStyle','-')
set(gca,'XGrid','off','YGrid','off','Fontsize',numsz)
xlabel('Years','Interpreter','LaTex','FontSize',labelsz)
ylabel('spread $\rho - r^*$ (\%)','Interpreter','LaTex','FontSize',labelsz)
xticks(tvec)
leg = legend('bad sunspot','good sunspot');
set(leg,'Interpreter','LaTex','Fontsize',legendsz)
legend boxoff
hold off
if savef == 1    
    paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
    paper_size = [3.0 2.0]; % [width height]
    paper_position = [0 0 3.0 2.0]; % [left bottom width height].
    paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
    fig = gcf;
    fig.PaperUnits = paper_unit;
    fig.PaperSize = paper_size;
    fig.PaperPosition = paper_position;
    fig.PaperOrientation = paper_orientation;
    cd(fdir)
    saveas(gcf,strcat('IRF_r_',cname,'_b',bval,'.pdf'))    
    cd(cdir)    
end

fig = figure(202202); clf
plot(tvec,100*IRF_nn_sB,'color',[0.8 0.4 0.4],'LineWidth',linesz,'LineStyle','--') ; hold on
plot(tvec,100*IRF_nn_sG,'color',[0.8 0.4 0.4],'LineWidth',linesz,'LineStyle','-')
set(gca,'XGrid','off','YGrid','off','Fontsize',numsz)
xlabel('Years','Interpreter','LaTex','FontSize',labelsz)
ylabel('issuance $n$ (\%)','Interpreter','LaTex','FontSize',labelsz)
xticks(tvec)
hold off
if savef == 1    
    paper_unit = 'inches'; % 'inches', 'normalized', 'centimeters', 'points'
    paper_size = [3.0 2.0]; % [width height]
    paper_position = [0 0 3.0 2.0]; % [left bottom width height].
    paper_orientation = 'portrait'; % 'portrait (default)', 'landscape'
    fig = gcf;
    fig.PaperUnits = paper_unit;
    fig.PaperSize = paper_size;
    fig.PaperPosition = paper_position;
    fig.PaperOrientation = paper_orientation;
    cd(fdir)
    saveas(gcf,strcat('IRF_n_',cname,'_b',bval,'.pdf'))    
    cd(cdir)    
end